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AN EFFICIENT NUMERICAL METHOD FOR SOLVING THE TIME-DEPENDENT 
COMPRESSIBLE NAVIER-STOKES EQUATI(»iS AT HIGH REYNOLDS NUMBER 


Robert H. MacComack 

Aaes Research Center, NASA 
Moffett Field. California, USA 


ABSTRACT 

This paper describes a new nuaerical nethod that has been used to drasti- 
cally reduce the computation time required to solve the Navier-Stokes equations 
at flight Reynolds nuid>ers. Though flows past complete aircraft configurations 
are still beyond our reach, the new method makes it possible and practical to 
calculate many important three-dimensional, high Reynolds number flow fields 
on today's computers. 


INTRODUCTION 

The Navier-Stokes equations adequately describe aerodynamic flows at 
standard atmospheric temperatures and pressures. If we could efficiently 
solve these equations, there would be no need for experimental tests to design 
fll^t vehicles or other aerodynamic devices. Unfortunately, analytic or 
closed-form solutions to these equations exist for only a few simple flow 
problems. During the past decade, the computer has been used to generate 
many new solutions. However, with existing numerical methods and computer 
resources, these solutions have been restricted to low Reynolds number or two- 
dimensional flows. This paper describes a new numerical nethod that has been 
used to drastically reduce the computation tine required to solve the Navier- 
Stokes equations at flight Reynolds numbers. Though flows past complete 
aircraft configurations are still beyond our reach, the new method makes it 
possible and practical to calculate many important three-dimensional, high 
Reynolds number flow fields on today's computers. 

The unsteady Navier-Stokes equations are difficult to solve because at 
high Reynolds numbers the magnitude of the inertial forces described by the 
hyperbolic terms of the equations are much larger than the viscous forces 
described by the parabolic terms. Conventional numerical methods, whether 
explicit or lsq)llclt, are severely restricted to small time steps by the 
stability or acctiracy requirements imposed by the hyperbolic terms; hence, 
many time steps may be required before the viscous effects can be determined. 
The present approach tine-splits the equations into a hyperbolic part and a 
parabolic part, solves the hyperbolic part by using a new explicit numerical 
method based on characterlsltics theory, and solves the parabolic part by using 
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« new efficient laplicit parabolic nethod. Both aethods are fully conservative 
and stable with tlM-steps orders of ae^ltude larger than those allowed by 
CFL (Courant. Friedrich, and Levy) stability criteria. 


DIFFEREMTIAl EQUATIONS 


The tiae-dependent coapressible Nevier'-’Stokes equations in two diaensions 
aay be written in conservation fora as 


3u.dF.dC_. 


( 1 ) 



in terw of density p, x and y velocity components u and v, viscosity 
coefficients X and u, total energy per unit voluae e, specific internal 
energy e, and coefficient of heat conductivity k. Finally, the pressure p 
is related to e and p by an equation of state, p(e,p), where 
e ^ (c/p) - [(u^ + v^)/2]. 


C(»fFUTATIOMAL MESH 

For calculating invlscid-vlscous interacting flows such as that sketched 
in Fig. 1 for a shock wave boundary-layer interaction, an efficient 



COTputatlonal aesh is used, shown in Fig. 2. It consists of a fine mesh near 
the wall for resolving the flow where viscous effects are important and a 
coarse mesh away from the wall where the flow is essentially inviscid. Addi- 
tional efficiency can be gained by stretching the fine mesh away from the wall 
to reduce the number of mesh points in regions requiring less resolution. 
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FDRMEK METHOD 

In 1970 a second-order accurate auaerical eethod was presented (1) 
for solving the unsteady conpressible Navler-Stofces equations for inviscid- 
viscoiis interaction probless. Since then, the method to be called "former" in 
this paper, has been significantly modified (^} to increase its accuracy 
and efficiency. A brief description of it follows. 

If the solution Ui j is known at time t ■= nAt at each mesh point 
(l,j), the solution at time t « (n + l}At is calculated by 




where ^(At) is a synmetric sequence of tiste-spllt, one-dimensional difference 
operators £j^(At^) and X^(Aty). For example. 






In this sequence the operator is called twice, each tine advancing the 

solution in time by At/2 by accounting only for the effect of the x-derivative 
in Eq. (1) on the solution. Similarly, the £y operator advances the solution 
by At once by accounting only for the effect of the y-derivatlve on the 
solution. The £y operator solves the time-split differential "equation" 


by first predicting a new value 

• "i.j 

and then correcting the predicted value. 


. 3g - 

from the current solution value U, , 

- 




+ u 


(p) _ ^ 

i.J Ay 




]) 


The corrected value then becomes the current value for the next split difference 
operator in the sequence. The operator is similarly defined. 
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The operators stable if 


At 

* |u| + c + (l/p){(2Y»i/PrAx) + t(-Ap)”2/Ay]} 

and 

At < ^ 

y |v| + c + (l/p){[(-Xy)“2/Ax] + (2Yu/PrAy)} 


idtere c is the speed of sound. Y is the ratio of specific heats of the gas, 
and Pr is the Prandtl nusd>er. 

For calculating an inviscid-viscous interacting flow on a t%io-mesh system 
shown in Fig. 2, typical operator sequences are (1) for all (i,j) in the coarse 
nesh 

%ihere At < min, ,{2 max At , max At }, and (2) for all (i,J) in the fine mesh 
“ i*j X y 

where m is the smallest Integer such that (At/m) < min. .(max At , 

2 max Aty}. 

For high Reynolds number calculations, the viscous region becomes very 
thin, requiring Ay near the wall to be very small. This causes Aty of the 
Xy operator also to be small and the integer m to be large. Values for m 
often exceed 100, requring a great amount of calculation time In the fine mesh. 

PRESENT METHOD 

Two new operators have been developed for replacing the Xy operator 
in the fine-mesh calculation 

The operator Xy_ solves for the inviscid (hyperbolic) terms of G. It 
is explicit, conservaflve, uses characteristic relations to predict convection 
and pressure fields, and is stable if 



Because the fine mesh covers only a thin layer of the flow adjacent to the wall, 
the normal to the wall velocity component v is very small relative to c; 
hence, the above stability bound for X_ is auch less restrictive than that 
for Xy. H 

Hie operator Xy solves for the viscous (parabolic) terms of G. It is 
laipllcit, conservative, requires no linearization, uses simple tridiagonal 
inversion instead of block trldlagonal inversion procedures, and is uncondi- 
tionally stable. 


The operator sequence for all (i,j) in the fine mesh for the present 
method is -im 






is 
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where now, beceuse of the greatly relaxed stability requirements, m is a s- II 
Integer usually equal to 2 in value. 


THE OPERATOR X 


The operator 



solves the split conservation lav "equation" 


dU . 


0 


where G„ 

a 



( 2 ) 


The convection velocity v and pressure p have been underlined for later 
reference. Conventional explicit difference methods fer this equation are 
limited to time steps 


At ^ 

jv| + c 


Violation of this condition, the CFL condition, is a violation of the domain of 
dependence principle and usually results in rapid numerical error growth. The 
operator Xy„ with time steps At < (Ay/|v|) circumvents the CFL condition 
by first preActing at time t -t- At the convection and pressure fields using 
characteristic relations which do not violate the domain of dependence prin- 
ciple and then by using the predicted values, values at time t, and the 
equations of motion to determine the time-averaged velocity and pressure 
gradients at each B»sh point (1,J) to be used to Integrate the above equation 
In a way similar to that described for the operator X^. 

To derive the required characteristic relations, we first write the above 
equation in nonconservation form 



If V is negligible compared to c, then the equations of motion become 


and 


If + YP - 0 


Using the relation c “ Vy we obtain th». following characteristic relations 
^ + pc 


<iC+ 


P 

dv n 4 . d 3 . 3 

^“0 where _ - ^ + c 


and 


dp dv - . d 3 3 

. PC ^ - 0 where " c 


These relations can be solved rapidly for large time steps on a character- 
istic set of mesh points such as that shown in Fig. 3. The characteristic mesh 
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• ORIGINAL MESH POINTS yj 

• CHARACTERISTIC MESH POINTS yjg 



Fig. 3 Characteristic laesh 


points {yj(.,Jc ■ 1,2,3 ...} are spaced according to the local speed of sound, 

although th'jy are shown equally spaced in the figure. The first characteristic 

nesh point is positioned c^tg/2 away from the wall, where c^ is the speed 

of sound at the wall and At^ - min^^ j(Ay/c). The jc + 1®^ mesh point is 

given by yjc+1 “ Tjc ®jc^^c '*ere c^^ is found by interpolation from 

values at the original set of mesh points {yj,j “ 1,2 The dependent 

variables Pj^, Vj^., and Pj^ are also obtained at each of the new mesh points 

by Interpolation. The new mesh has the following property: a sound signal 

traveling either to the right or left can leave any point in the mesh at time 
t and arrive at its next neighbor at t + At^., and at its next neighbor at 
t + 2Atc, etc. This property can be used to solve quickly for v and p at 
time t* - t + kAtg, as shown in Fig. 4, where k is the largest integer such 
that kAt < mlnj .(Ay/|v|). Approximating the characteristic relations by 

fjc - pjo-k * ‘’jc'jo*'';. - ’jc-k’ ■ “ 
and 

**jc *’jc+k *’jc‘^jc^^jc '’^jc+k^ 

wf obtain 


* 

Plc-k ^ ‘'lc+k 

+ D , C . 

’jc 

" 2 

*^jc jc 

* 

'ic-k * ’jekk 

+ !ic-k ■ 

jc 

2 



^1c-k ~ ^1c4k 


- P. 


Pot characteristic paths reflecting from the wall, jc - k i 0, and “d 

are replaced by Pj_(jc-k) ~'^l-(jc-k) ' '^*“P®ctlvely. 
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An altemrtlve solution technique and the one used by the present anthod 
is to replace the point values by integrated space-averaged values. Thus, 


and 


where 


* ^-1c ^+1c . 

Pjc“"2— Vjc 


-ZlS: ±i£. 

2 


* . ±l£ + _li£ +lc 


P --i— fl 

-jc 2k + 1 Pjc-A • *^+jc 2k + 1 ^ 

2k 2k 

-jc " 2kTT ^ '^jc-i • '^+jc “ 2T+T 2 


jc+i ’ 


Jc+i 


Once the Integrated averages have been detemined for point jr, only four 
additions and four subtractions are required to obtain the values for point 
jc + 1 for any k. 


Using the values of v and p at times t and t* and the equations of 
motion, we can determined the time-averaged velocity and pressure gradients 
responsible for the change In the solution 


and 


where At ■ kdt ■ t 
c 
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Th«M an the velocity and pnasun gndlanta which will ha uaad to Integrate 
I4. (2). Ue fltat nuaerlAally Integra the grwllanta la y to determine v 
nd p at aach of the original eat of Mah polnta {y^tJ * l>2t3 ...} 


and 


vf 




p^-p. 


& 


0 'i' ^ 


where p^ ir preaaure at the mil and la found by a characteristic relation 
aiallar to the ones given earlier. Equation (2) can now be Integrated in time 
by At ■ kAtc in the fine mesh in the sme manner as described for the former 
method except that the convection velocity and pressure (the underlined terms 
of Gg) are nplaced by their time-averaged values dcrarmined above in both the 
predictor and corrector steps. 

raE (AERATOR £„ 


The operator £y solves the split conservation law "equation" 


0 where 6_ 




(3) 


This equation is parabolic and from domain of dependence considerations it is 
most appropriately solved using an implicit numerical method. Examples of such 
methods applied to the following model parabolic equation 

^ + [small terms] 

are (a) Crank-Nicolson 


0+1 « “ j. / n , n . n \ 
u . “ U 4 + — I u . , - 2u. + u . . , 1 
j J 2Ay2\3“l i 

( n+l , n +1 n+l\ 


2 Ay2 

+ 


2Ay2 

■f At (small terms] 


and (b) Laasonen 


u^^ • u“ + ^ ^u^l - 2u^^ + Uj+lj + [small terms]© 

The above difference equations can be efficiently solved by a simple tridlago- 
nal inversion. 

Equation (3) can be written 

ll-* 

|£ . i 4 i^iMC9y/ax)i| 
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Th« flrit of the ebove equations, the continuity equu^n without convec*- 
tion, indlcstes that density is stationary in time during ^m.s time-split 
interval and is thus trivially solved. The second uid third equations, the x 
and y aoaentua equations, have essentially the same form as the model parabolic 
equation and can be efficiently solved using two simple trldlagonal inversions. 
Unfortunately, the fourth equation, the total energy equation, does not have 
the same form as the model equation because th^ solution variable 
e 4- [(u^ 4- v^)/2] does not appear on the right-hand side in a second derivative i| 
with respect to the y term, although parts of it do. To avoid using a more 
costly block trldlagonal procedure for solution, we split the energy equation 
into three equations, a u^ kinetic energy equation, a v^ kinetic energy 
equation, and an Eternal energy equation. We obtain the kinetic energy 
equations by multiplying the x and y snmenCum equations by u and v, 
respectively. The internal energy equation is obtained by subtracting the two 
kinetic energy equations from the total energy equation. The resulting equa- 
tions in model parabolic form are 


^ . i 9[M(3u^/9y)] ^ II 3[2mu(3v/3xJ] _ ^ 
dt '' oy \p 5y 1 



The u^ and v^ kinetic energy equations form the same trldlagonal coefficient 
matrices as the u and v momentum equations; hence, solution of the three 
energy equations requires the inversion of only one additional tridiagonal 
matrix. 


The operator Xy solves Eq. (3) using either the Crank-Micolson or the 
Laasonea method. The^second-order accurate Crank-Nieolson method is preferred 
unless the coefficient pdt/dy^ is large enough that the numerical amplifica- 
tion factor of the method approaches -1, at \diich point the Crank-Nlcolson 
method is known to believe erratically. For large values of y&t/Ey^, the first- 
order accurate Laasonen method is used. The variations of p, M, and A in y 
cause no difficulty for either method. The variations of u, >■, and the 
[small terns] on the right cn be accounted for by solving the equations twice, 
first with the present values of y, X, and [small terms] and then with the 
newly calculated values. Averaging the two solutions time-centers the differ- 
ence equations. The assumption that the mixed derivative terms can be treated 
as small terns has caus^ no numerical difficulties. The present method first 
solves the x and y momentum equations, uses the present and new u and v 
solutions to evaluate the dissipation terms T^ and Tj (not always small) in 
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«ith«r a Crank-Nlcolsoii or Laasonon oannar, and than aolvas tha v^, rad c 
raargy aquatlona. Errora aada In tha and t aquations bacausa of 

larga T< rad T> cracal axactly when tha aol,<tlona ara coa^lnad to obtain 
a ■ p{e + ((u^ + v^)/2]}. 

Tha oparatoT ty^ invarta thraa aiapla tridlagonal confidant siatricas 
and aolvas flva ayataas of aquations. Using tha convantion of counting only 
nultlplicationa and dlvialraa (3), tha prsasnt nathod raquicas 21M • 16 arlth- 
Mtic operations » whare N is tha order of the tridiagonal systca. Conven- 
tional adhods for solving Eq. (3) invert one block tridlagonal aatrlx ot 
order ’Ipwith block eleaent utricej of order thraa and require 108N - 18 
operations if the Inverses of tha diagonal block elenent oatricas are computed, 
or approxiaately 4SH operations if tha Inverses ara not computed. 

RESULTS 

Several shock wave boundary-layer interaction problraw vara calculated 
using the foraar aathod {2) present method. For each calculation the 

flow was at Madi 2, and a shock wave Incident upon a flat plate increased tha 
pressure by a factor of 1.4. Molecular viscosity was calculated using 
Sutherland's formula, and turbulent eddy v^uoslty was calculated using a 
Cebaci-^l^th turbulence model. Figure 5 cffl^ares Che results of both methods, 
experismnt and boundary-layer theory (Ref. using Crocco's method) for 
a separated laminar botmdary layer at a ReynoldB*~ninAar of 2.96 x 10^. The 
results of the two amthods agree well everywhere with the exception of pressure 
at the leading edge tfig* 5(a)] which will be discussed herein. The present 
method required 1.0 min of CDC 7600 computer tliM, and the former method 
required 6.4 min. Figure 6 compares the results of both methods for a turbu- 
lent boundary-layer Interaction at a Reynolds number of 2.96 x 10^. Again 
Che results of the two methods agree well except at Che leading edge. Fig- 
ure 6(c) displays the pressure profiles at the leading edge predicted by both 
methods. The present method is far more sensitive to the presence of the 
leading edge singularity chan the former method and appears to be correctly 
approaching the theoretical value of pressure behind a Mach 2 normal shock 
nave. The present method required 1.4 min of computation time, and Che former 
method required 58. S min. Figure 7 compares Che computing times of the former 
and present methods for a wide range of Reynolds numbers. 
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Fig. 5 Continued 
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(c) Leadlng-ed^ singularity 
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(d) ValQCity profile ahead of interaction 


Fig. 6 Continued 
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(e) Velocity profile Interaction region 


y/HF 



(f) Velocity profile aft of interaction 
Fig. 6 Concluded 
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Fig. 7 CoBparlson of foraer and present method computing times on the 

COC 7600 vs. Reynolds ntmber for several shock boundary-layer interaction 
calculations 


COHCLUSION 

The present ^thod has reduced the computation time by one and two orders 
of magnitude from that required previously to solve for the interaction of a 
shock wave with a boundary layer on a flat plate. The method has been applied 
with similar success to several other invisc id-viscous Interaction problems, 
including flows past con|>resslon comers, wavy walls, axisymmetric channels, 
and blunt-nosed lifting airfoils, at Mach nui^ers as high as 8.5, with shock 
waves strong enough to increase surface pressure by a factor of 80, and at 
Reynolds numbers as high as 10^. 
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